%--------------------------------------------------------------------------
%   天线方向图函数
%   qwe14789cn@gmail.com
%   20180402
%--------------------------------------------------------------------------
%   坐标系描述
%   阵列法线方向指向x轴正向，阵元分布阵列面按照y轴分布
%   o 表示阵元
%   ->x 表示法线指向
%   ->  表示子阵元天线方向图主轴指向
%   各个区间表示目标指向角度
%--------------------------------------------------------------------------
%
%                               | y
%                               |
%          [-90 -180]           |        [0 -90]
%                               o  ->
%                               o  ->
%                               o  ->
%   ------------------------------------------------------------> x
%                               o  ->
%                               o  ->
%                               o  ->
%          [90   179]           |        [0  90]
%                               |
%                               |
%
%--------------------------------------------------------------------------
%   输入:
%       所有阵元的空间坐标，格式如下
%       array = |x1 x2 x3 ...|
%               |y1 y2 y3 ...|
%               |z1 z2 z3 ...|
%       lambda      射频信号波长              单位 米
%       w           导向矢量+窗函数，按照阵元的顺序写，角度和幅度 复数
%       sub_angle   对应array中的子阵元主轴中心指向角度，例如圆阵指向角度
%                   并不像线阵一样都是0度,而是如[-4 -2 0 2 4]这样
%       sub_pattern  子阵元天线方向图,-180~179度,间距1度,第二列为振幅增益
%       sub_pattern =    |-180 A1|
%                       |-179 A2|
%                       |-178 A3|
%                       | ... ..|
%                       | 179 An|
%   输出：
%       theta   各个角度坐标
%       E       矢量电场强度
%--------------------------------------------------------------------------
%   1.线阵4阵元无导向矢量
%   输入
%       目标距离
%       输入阵列形态
%       射频波长
%--------------------------------------------------------------------------
%   example 1
%--------------------------------------------------------------------------
% shape = [0 0 0 0;-1.5 -0.5 0.5 1.5;0 0 0 0];
% lambda = 2;
% [theta,E] = array_pattern(shape,lambda);
% figure(1)
% plot(theta,pow2db(abs(E).^2));grid on;
%--------------------------------------------------------------------------
%   2.线阵5阵元+导向矢量
%   输入
%       目标距离
%       输入阵列形态
%       射频波长
%       导向矢量角度
%       导向矢量
%--------------------------------------------------------------------------
%   example 2
%--------------------------------------------------------------------------
% shape = [0 0 0 0 0 0;-2.5 -1.5 -0.5 0.5 1.5 2.5;0 0 0 0 0 0];
% lambda = 2;
% w_angle = -20;
% w = exp(-1i*2*pi/lambda*1*sind(w_angle)).^(0:5)';
% [theta,E] = array_pattern(shape,lambda,w);
% figure(2)
% plot(theta,pow2db(abs(E).^2));grid on;
%--------------------------------------------------------------------------
%   3.多阵元+导向矢量+窗函数
%   输入
%       目标距离
%       输入阵列形态
%       射频波长
%       导向矢量角度
%       导向矢量
%       子阵元指向角度
%       子阵元天线方向图
%--------------------------------------------------------------------------
%   example 3
%   暂无例子
%--------------------------------------------------------------------------
function [theta_angle,E] = array_pattern(array,...
                                        lambda,...
                                        w,...
                                        sub_angle,...
                                        sub_pat)
c = 299792458;                                                              %光速
k = 2*pi/lambda;                                                            %波数
step = 0.1;
theta_angle = -180:step:179;                                                     %输出角度
tgt_range = 10000*lambda;
%--------------------------------------------------------------------------
%   输入变量数量判断
%--------------------------------------------------------------------------
if nargin <= 2
    w = ones(1,length(array));
    sub_angle = zeros(1,length(array));
    sub_pat =[-180:179;ones(1,360)]';
elseif nargin <=3
    sub_angle = zeros(1,length(array));
    sub_pat =[-180:179;ones(1,360)]';
end
%--------------------------------------------------------------------------
%   子阵天线方向图插值到0.1
%--------------------------------------------------------------------------

sp(:,1) = theta_angle;
sp(:,2) = interp1(sub_pat(:,1),sub_pat(:,2),theta_angle);
point = 1;
for degree = -180:step:179                                                  %目标电磁波角度穷举
    tgt = [tgt_range*cosd(degree);tgt_range*sind(-degree);0];               %目标空间坐标
    for index = 1:length(array(1,:))                                        %子阵元循环
        A_sub_angle = circshift(sp,round(sub_angle(index)/step));
        E(point,index) = exp(-1j*k*norm(tgt - array(:,index))).*...
                        w(index).*...
                        A_sub_angle(point,2);
    end
    point = point + 1;
end
